# Make graphs for figure 1

# Load Models and Store AUCs
countries <- c('indo', 'colo')
dvs <- c("any", "high", "spike")
titles <- c(any = "Any violent event",
            high = "Five or more violent events",
            spike = "1SD Increase in events")

for (country in countries) {
  for (v in dvs) {
    source(paste("config_code/",country,"_data_setup.R",sep = ""))
    source(paste("config_code/",country,"_predictor_vars.R",sep = ""))
    models <- names(figures_1_2)[-1]

    # Get population benchmark dataset
    pop_auc <- array(NA,
                           dim = T)
    filename <- paste(modeldir,"/",
                      country,
                      "_ebma_",
                      v,
                      "_population_cross.RData",
                      sep = "")
    load(filename)
    predictions <- as.vector(ebma.results[[1]]$fit.oos)
    realizations <- as.vector(ebma.results[[1]]$actual.oos)
    pop_auc <- roc(response = realizations,
                   predictor = predictions)[['auc']]
    
    auc_array <- array(NA,
                       dim = c(length(models)),
                       dimnames = list(models))
    for (x in models) {
      filename <- paste(modeldir,"/",
                        country,
                        "_ebma_",
                        v,
                        "_",x,"_cross.RData",
                        sep = "")
      load(filename)
      predictions <- as.vector(ebma.results[[1]]$fit.oos)
      realizations <- as.vector(ebma.results[[1]]$actual.oos)
      auc_array[x] <- (roc(response = realizations,
                              predictor = predictions)[['auc']]
                          - pop_auc)
    }

    # Make Graph
    var.order1 <- order(auc_array[names(figures_1_2_fixed)[-1]],
                        decreasing = FALSE)
    var.order2 <- order(auc_array[names(figures_1_2_slow)],
                        decreasing = FALSE)
    var.order3 <- order(auc_array[names(figures_1_2_annual)],
                        decreasing = FALSE)
    ranks1 = rep(NA,length(var.order1))
    ranks1[var.order1] <- ((1:length(var.order1))
                           +length(var.order2)
                           +length(var.order3)
                           +ifelse(length(var.order2)>0,
                                   2,1))
    ranks2 = rep(NA,length(var.order2))
    ranks2[var.order2] <- 1:length(var.order2)+length(var.order3)+1
    ranks3 = rep(NA,length(var.order3))
    ranks3[var.order3] <- 1:length(var.order3)
    ranks <- c(ranks1,ranks2,ranks3)
    varsets.total <- c(names(figures_1_2_fixed)[-1],
                       names(figures_1_2_slow),
                       names(figures_1_2_annual))
    pdf(paste("graphs/",
              country,
              "_AUCs_",
              v,
              "_ebma_popplus_cross.pdf",
              sep=""),
        width = 7, 
        height = 7)
    par(mar=c(5, 12, 2, 1) + 0.1)
    plot(auc_array[varsets.total],
         ranks,
         xlim = c(-.1,.1),
         ylim = c(0,max(ranks)+ifelse(length(var.order2)>0,
                                      2,1)),
         ylab = "",
         xlab = "Performance improvement over population (AUC)",
         yaxt='n',
         main = "")
    for (i in 1:length(ranks)) {
      lines(x = c(-1,1),
            y = c(ranks[i],ranks[i]),
            col = "gray95")
    }
    points(auc_array[varsets.total],
           ranks)
    axis(2,
         at = ranks,
         labels = labels[varsets.total],
         las = 2)
    axis(2,
         at = c(length(ranks)+ifelse(length(var.order2)>0,
                                     3,2),
                length(ranks3)+1,
                length(ranks2)+length(ranks3)+ifelse(length(ranks2)>0,2,1)
                ),
         labels = c("Time-Invariant",
                    "Time varying (Annual)",
                    ifelse(length(ranks2)>0,"Slow-Moving","")),
         las = 2,
         font=4 )
    lines(x = c(0,0),
          y = c(-1,length(ranks)+ifelse(length(ranks2)>0,5,3)),
          lty = 3)
    dev.off()
  }
}
